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ABSTRACT 


The  effects  of  the  rotation  stress  mechanism  and  the  penetration  of  short- 
wave radiation  below  the  sea-surface  are  examined  in  determining  a  one- 
dimensional  equilibrium  mixed  layer  depth.  Starting  with  the  Obukhov-scale 
equilibrium  theory  for  the  surface  ocean  boundary  layer,  a  revised 
equilibrium  theory,  which  includes  rotation  stress  and  radiation  effects,  is 
presented.  This  new  theory  is  applied  using  climatological  boundary 
conditions  for  the  tropical  Atlantic,  and  the  results  are  compared  with  the 
observed  climatological  mixed  layer  depth. 

In  general,  the  response  of  the  revised  model  is  an  improvement  over  the 
Obukhov  theory  alone.  Because  the  quality  of  the  results  are  limited  by 
uncertainties  in  the  boundary  conditions,  no  detailed  evaluation  of  the  model 
response  is  justifiable.  However,  it  is  concluded  here  that  the  physical 
mechanisms  of  rotation  stress  and  penetration  of  radiation  are  important  in 
determining  a  steady-state  equilibrium  depth  of  mixing  for  the  tropical 
Atlantic. 
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.    INTRODUCTION 


The  study  of  the  oceanic  turbulent  boundary  layer  is  a  relatively  recent 
field  in  Physical  Oceanography.  Since  the  pioneering  work  of  Kraus  and 
Turner  (1967),  a  significant  effort  has  been  made  in  understanding  and 
modelling  the  dynamics  of  turbulent  mixing  in  the  oceans.  However,  the  basic 
contributions  made  by  those  authors  still  provide  the  theoretical  foundation 
for  much  of  the  current  research. 

The  oceanic  boundary  layer  or  mixed  layer  is  considered  in  this  study  to 
be  the  fully  turbulent  region  that  is  bounded  above  by  the  air-ocean 
interface,  and  it  is  where  the  temperature  and  salinity  are  usually  observed 
to  be  fairly  well  mixed.  Below,  the  mixed  layer  is  assumed  to  be  bounded  by  a 
dynamically  stable  thermocline. 

The  study  of  the  top  few  tens  of  meters  of  the  ocean  is  of  considerable 
scientific  interest.  Ekman  pumping  effects,  originating  in  the  mixed  layer  by 
the  action  of  wind  forcing,  considerably  influences  the  dynamics  of  the  lower 
levels  and  the  interior  circulation  in  general.  Also,  the  interaction  between 
the  oceanic  and  atmospheric  boundary  layers  is  an  essential  mechanism 
which  must  be  considered  when  making  medium  and  long-range  weather 
forecasts,  since  a  large  part  of  the  atmospheric  energy  supply  comes  from 
the  heat  exchanged  between  these  layers  .  Additionally, this  is  the  region  of 
primary  biological  productivity,  which  is  of  significant  economic  importance. 
An  important  military  application  is,  on  the  other  hand,  the  modelling  of 
accoustic  propagation  in  the  oceans. 
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Most  of  the  physics  behind  existing  models  of  the  oceanic  mixed  layer  are 
based  on  the  flux  form  of  the  Navler-Stokes  equations  of  motion,  with  the 
Bousstnesq  approximation.  One-d1mens1onal  versions,  like  the  one  which  Is 
treated  In  this  study,  further  assume  the  horizontal  gradients  of  the  mean 
fields  to  be  negligible.  For  a  steady-state  situation,  this  will  lead  to  the 
basic  energy  balance,  as  stated  by  Kraus  and  Turner  (1967),  between  the 
work  done  by  the  wind  stress  and  the  surface  heating,  in  defining  an 
equilibrium  depth  of  mixing.  However,  a  usually  neglected  term  in  the 
turbulent  kinetic  energy  budgets,  the  planetary  rotation  term,  was  recently 
examined  by  Garwood  eta/.  (1985a).  They  suggest  a  new  formulation  for  the 
equilibrium  mixed  layer  depth,  which  basically  describes  the  interaction 
between  the  zonal  Reynolds  stresses  and  the  northward  component  of  the 
planetary  rotation.  This  new  formulation  might  partially  explain  the  zona! 
variation  of  the  observed  mixed  layer  depth  in  the  Tropical  Pacific  and 
Atlantic  Oceans. 

The  goal  of  this  study  is  to  apply  the  Kraus  and  Turner  (1967)  and  the 
Garwood  et  al.  (1985b)  one-dimensional,  steady-state  models  on  a  relatively 
broad  grid  with  realistic  boundary  conditions  in  the  Equatorial  and  Tropical 
Atlantic  Ocean,  and  to  compare  the  resulting  diagnostic  mixed  layer  depths 
with  the  observed  mixed  layer  depths.  In  this  way,  insight  may  be  gained  Into 
the  relative  importance  of  the  physical  mechanisms  involved.  Also,  we 
should  be  able  to  Identify  the  regions  and  time  of  the  year  where  the  turbulent 
boundary  layer  can  be  represented  realistically  by  such  simple  steady-state 
equilibrium  models. 


II.     THE   DYNAMICS   OF   MIXING 


A.  THE  TURBULENT  KINETIC  ENERGY  BUDGET 

In  the  presence  of  a  constant  wind  stress  and  a  downward  surface  heat 
flux,  the  steady-state,  one-dimentional  turbulent  kinetic  energy  (TKE) 
budget  is,  assuming  no  vertical  advection  : 

0=d/dt  [E/2]= 


-u'w"  du/dz  -v'w'  dv/dz  +  b'w'-d/dz  [w'(E/2)+  w'p'/p]  -  d,  (2.1) 

where  E/2  is  the  TKE,  -u'w'  du/dz  -v'w'  dv/dz  represents  the  production  of 
kinetic  energy  associated  with  the  vertical  shear  in  the  horizontal  flow 
induced  by  the  wind  stress,  b'w'  is  the  buoyancy  flux,  -d/dz  [w'(E/2)+w'p'/p] 
is  the  turbulent  diffusion  and  d  is  the  dissipation.  With  no  loss  of  generality, 
we  shall  assume  the  x-axis  to  be  oriented  downstream,  sothatv=0. 
1 .  Shear  Production  of  Turbulence 

Let  u(z)  be  some  monotonic  function  of  z  between  the  surface  ,  z=0, 
and  the  bottom  of  the  mixed  layer,  z=-h.  For  the  sake  of  simplicity,  it  is 
further  assumed  that  u(-h)=0,  since  no  entrainement  is  allowed. 

Consider  now  (Fig.  la)  that  some  particle  initially  at  1  is  displaced  to 
some  position  2  by  means  of  a  positive  turbulent  vertical  velocity  w'. 
Because  du/dziO,  the  particle  will  acquire  a  negative  turbulent  horizontal 
velocity    u'    and,  then,     u'w'    will  be  less  tha     zero  at  2.  Alternatively,   a 
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particle  Initially  at  1 ,  which  1s  displaced  to  3  by  means  of  a  negative  vertical 
turbulent  velocity  w',  will  acquire  a  positive  horizontal  turbulent  velocity  u" 
and,  again,  u'w'  will  be  negative.  In  both  situations,  the  shear  production 
term  -u"w"  du/dz  will  be  positive,  and  when  this  physical  process  is 
averaged  over  the  whole  turbulent  layer,  it  will  tend  to  increase  the  TKE  at 
the  expense  of  the  mean  kinetic  energy.  Formally, 

-uV  du/az  >  0. 

Evaluated  over  the  mixed  layer,  the  shear  production  will  be: 

o 

G  h  (  -u'w'  du/dz  dz  =  <-u'w'>  [u(0)-u(-h)],  (2.2) 

-h 

where  <  >  represents  a  vertical  average.  At  the  bottom  of  the  mixed  layer, 
and  under  the  assumption  of  no  entrainement,  -iTw'(-h)  is  zero.  Then, 
disregarding  nondlmensional  constants,  <-i?w'>  Is  of  order  -i7w'(0).  Also, 
we  have  already  assumed  u(-h)  to  be  zero.  Thus, 

G  =  -CTw'CO)  u(0).  (2.3) 

Furthermore,  the  Reynold  stress  -i7w'(0)  can  be  written 

|  -lTw'CO))  =  |T|/p  ,  (2.4) 
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where  p  Is  a  typical  sea  water  density  and    i  Is  the  surface  wind  stress. 
Finally,  (2.3)  Is  rewritten  as 


G  =   u(0)|t|/p  . 


(2.5) 


2.  The  Buovancv  Flux 

Under  the  Influence  of  surface  heating,  there  will  be  a  tendency  for 
stratification,  and  mixing  will  only  occur  in  the  presence  of  wind  stirring.  In 
this  case,   the  less  dense  particles   near  the   surface   will    mix  with  the 
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Fig.  1.  (a)  Shear  production  of  turbulent  kinetic  energy,  In  the  presence  of 
wind  stress,  (b)  Buoyant  damping  of  turbulent  kinetic  energy,  in  the 
presence  of  surface  heating. 
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underlying  water  only  at  the  expense  of  the  wind-generated  turbulence,  thus 

converting  TKE  Into  potential   energy  (Fig.    lb).   This   means  that,   on  the 

average,  the  buoyancy  flux   -b'w'   will  be  a  sink  term  In  the  TKE  budget.  The 

equilibrium  between  the  effects  of  wind  mixing  and  surface  heating  1s  Indeed 

the  fundamental  energy  balance  for  any  steady-state,  one-dlmentlonal  mixed 

layer  model,  as  long  as  advectlon  Is  neglected. 

To  relate  the  buoyancy  flux  to  the  surface  heat  flux   Q0    we  apply  the 

linearized  equation  of  state  and  obtain: 

b'w'  =  a  g  TV,  (2.6) 

where  a  is  the  thermal  expansion  coefficient,  g  1s  the  acceleration  of 
gravity  and  TV  1s  the  temperature  flux.  It  can  further  be  shown  that, 
assuming  vertical  homogeneity  in  the  mean  temperature  T,  TV  will  be  linear 
over  the  mixed  layer,  and 

TV(0)--Qo/(  pCp),  (2.7) 

where  Q0,  the  surface  heat  flux  is  positive  downward. 

Integrating  the  buoyancy  flux  over  the  mixed  layer, 

o o  

BH  J  b'w'  dz  =  ag  J  TV  dz  =  agh  [TV(0)+TV(-h)]/2,  (2.8) 

where  h  is  the  mixed  layer  depth.  Because  no  entrainment  Is  assumed, 
TV(-h)=0.  Thus, 
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B  =   -aghQ0/(2pCp).  (2.9) 

This  simple  model  for  the  vertically  Integrated  buoyancy  flux  also 
requires  the  assumption  that  the  net  solar  radiation  be  completly  absorbed  at 
the  surface.  We  will  see  in  the  following  sections  that  this  assumption  is 
incorrect  when  the  vertical  scale  of  penetration  of  radiation  below  the 
surface  is  comparable  to  the  vertical  scale  of  turbulent  mixing.  The  effect  of 
the  penetration  of  radiation  is  expected  to  be  relevant  in  optically  clear 
oceanic  waters  and  in  summertime  conditions,  when  the  mixed  layer  is 
relatively  shallow.  In  those  situations,  the  net  downward  buoyancy  flux  may 
be  significantly  reduced,  causing  the  resultant  depth  of  mixing  to  be 
significantly  increased. 

3.  Turbulent  Diffusion 

Although  the  diffusion  of  turbulence  is  an  essential  mechanism  in 
transporting  TKE  downwards  from  the  upper  levels,  it  does  not  constitute,  by 
itself,  a  source  or  sink  of  turbulence,  which  means  that  its  net  value  over 
the  mixed  layer  must  vanish.  Formally, 

(  -d/dz  lw'(E/2)+Wp7p]  dz  =  0.  (2.10) 

4.  Dissipation 

Early  mixed  layer  models,  including  the  Kraus  and  Turner  (1967) 
model,  recognized  the  existence  of  dissipation.  However,  they  considered  its 
influence  to  be  negligible,  given  the  relatively  large  length  scale  of  turbulent 
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mixing.  Indeed,  a  less  careful  scale  analysis  of  the  TKE  budget  seems  to 
support  the  approximation.  More  recent  models,  however,  include 
dissipation  In  the  mechanical  energy  budget,  either  as  a  fraction  of  the  wind- 
-stress  production  or  as  a  fraction  of  the  net  TKE  production.  Although 
viscous  dissipation  may  be  neglected  when  dealing  with  turbulent  eddies  of 
length  scales  of  order  of  the  mixed  layer  depth,  the  dynamic  interaction 
between  those  eddies  produces  other  eddies  of  smaller  and  smaller 
dimensions,  where  the  role  of  dissipation  may  finally  dominate.  This  fact  is 
easily  understood  with  the  analysis  of  a  spectrum  of  turbulence,  where  it  can 
be  verified  that  the  mechanical  energy  drops  rapidly  above  a  limiting 
frequency  (Tennekes  and  Lumley,  1962). 

If  such  a  dissipation  mechanism  is  included  in  the  TKE  budget  by 
parameterization,  for  example  as  a  fraction  of  the  TKE  production,  we  will 
see  in  the  following  sections  that  the  Kraus  and  Turner  model  is  still  usable, 
since  no  adltional  degree  of  freedom  is  required.  That  was  the  approach  of 
Gelsler  and  Kraus  (1969),  Niller  (1975)  and  others. 

For  now,  we  will  just  define  the  vertically  integrated  dissipation  D  as: 


o 

D^Jddz.  (2.11) 


-h 


B.   EQUILIBRIUM  DEPTH  OF  MIXING 

From  the  results  of  the  previous  section,  the  vertically  integrated  TKE 
budget  can  be  written  as 
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B  +  G-D  =  0, 


(2.12) 


where 
and 


B  =  -  agh  Q0/(pCp) 
6  =  u(0)|t|/p  . 


Work  by 
wind  stress 


V 


Mean  Kinetic 
Energy 


Shear 
production 


\S 


Surface 
heat  flux 


Potential 
Energy 


/\ 


Buoyant 
damping 


Horiz.  Turbul 
Kinetic  energy 


Redistribution 

i         :> 


c 


Yert   Turbul 
Kinetic  energy 


^O^  Dissipation  // 

(  Heat  ) 


Fig.  2.  The  one-dimensional  mixed  layer  equilibrium.  The  "redistribution* 
effects  include  the  pressure  redistribution  and  the  rotation  stress 
mechanism. 
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Equation  (2.12)  illustrates  the  basic  steady-state  balance  between  the  work 
done  by  the  wind  stress  and  the  buoyant  damping  (Fig.  2).  Note  that,  because 
of  turbulent  mixing,  the  mean  flow  will  tend  to  be  vertically  homogeneous, 
except  near  the  surface,  where  turbulence  Is  generated,  and  near  the  bottom 
of  the  mixed  layer,  where  it  falls  to  zero.  Also,  the  mean  temperature 
profile  will  tend  to  be  homogeneous.  Thus,  the  mixed  layer  potential  energy 
will  approach  its  theoretical  maximum  value. 

In  Fig.  2  an  aditional  mechanism  is  Illustrated.  This  processs,  called 
"redistribution",  does  not  appear  in  the  TKE  budget,  since  its  constituent 
terms  sum  to  zero  when  the  three  spatial  components  of  the  kinetic  energy 
budget  are  added.  Redistribution  consist  of  two  parts:  the  pressure 
redistribution  and  the  planetary  rotation.  Their  roles  are  to  exchange  TKE 
among  the  x-y-z  directions,  tending  to  restore  the  system's  Isotropy.  In 
section  III,  the  planetary  rotation  term  will  be  shown  to  have  a  fundamental 
role  in  determining  the  equilibrium  mixed  layer  depth  for  the  Garwood  et  al. 
(1985a)  steady-state  model. 

1.  TheObukhov  Length  Scale 

We  have  seen  that,  if  dissipation  is  parameterized  in  terms  of  the 

other  basic  turbulent  variables,  the  steadiness  of  the  mixed  layer  depth  will 

depend  on  the  surface  boundary  conditions,  i.e. ,  the  surface  wind  stress  and 

the  surface  heat  flux.  This  is  equivalent  to  setting  G/B  =  constant. 

We  will  now  define  a  velocity  scale   u#  ,  which  is  representative  of  the 

typical  turbulent  speed  of  the  largest  eddies  in  a  wind-driven  mixed  layer: 

u»2.|T|/p.  (2.13) 
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Because  u(0)    depends  linearly  on  the  wind  s  eed  and,  thus,  on  the  surface 
wind  stress,  the  shear  production  term  G  In  equation  (2.5)  can  be  scaled  as 

G~u*3.  (2.14) 

Combining  (2.9)  and  (2.12-14),   and  setting     -G/B  =1,  the  Obukhov 
length  scale  is  defined  as 

L  =  h  =  2pCpu*3/  (agQ0).  (2.15) 

Physically,    L    represents  the  scale  of  maximum  depth  of  turbulent  mixing, 
given  the  surface  boundary  conditions  of  wind  stress  and  heat  flux,  for  an 
equilibrium  mixed  layer  with  no  advection  present. 
2.  The  Kraus  and  Turner  Steadv-State  Model 

Based  on  equation  (2.12),  Kraus  and  Turner  (1967)   considered  the 
dissipation  to  be  negligible,  so  that 

B  +  G  =  0.  (2.16) 

Also,  the  vertically  integrated  shear  production  G  can  be  written 

G  =  u*3=  (N/p)3/2  (2.17) 

or  G=  (paCDu2/p)3/2  f  (2.18) 
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where   C&   is  the  drag  coefficient,    pa    is  the  air  density  and   U    is  the  wind 

speed  measured  at  10  meters  above  the  surface.  Substituting  (2.9)  and 
(2. 18)  into  (2. 16),  the  equilibrium  depth  of  mixing  is  written 

h  =  Ci  (  p  Cp/  a  g  Q0  )(  Pa  CD  u2/p)3/2  ,  (2.19) 

where  C;   is  a  dimensionless  constant  of  order  1 . 

More  recent  versions  of  the  Kraus  and  Turner  model  assume  the 
dissipation  D  to  be  proportional  to  the  total  TKE  production  B+G,  so  that 
equation  (2. 12)  is  written  as 

(K+  1)(  B  +  G  )  =  0,  (2.20) 

where  K  is  a  dimensionless  constant.  As  stated  earlier,  this  will  lead  to  the 
same  basic  result  of  equation  (2. 19),  since  no  aditional  degree  of  freedom  is 
included. 
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III.     THE   REVISED   STEADY-STATE   MODEL 


As  noted  in  section  I,  the  dynamics  of  the  oceanic  mixed  layer  in  the 
equatorial  and  tropical  regions  are  not  fully  understood.  One  of  the  most 
apparent  but  least  understood  features  is  the  zonal  variation  of  the  mixed 
layer  depth. 

Some  of  the  discrepancies  between  the  observed  mixed  layer  depths  and 
the  predictions  made  by  existing  steady-state  and  one-dimensional  models 
might  be  explained  by  the  exceptionally  strong  vertical  and  horizontal 
advection  which  are  characteristic  of  the  tropical  oceans.  Also,  the  effect  of 
penetration  of  short  wave  radiation  below  the  surface  is  believed  to  have 
some  importance  for  relatively  shallow  mixed  layers,  especially  in  the 
presence  of  large  net  solar  irradiance.  Another  physical  mechanism, 
examined  by  Garwood  et  ah  (1985a),  involves  the  interaction  between  the 
zonal  surface  wind  stress  and  the  planetary  rotation,  and  it  may  partially 
explain  the  zonal  dependence  of  the  steady-state  mixed  layer  depth. 

It  is  beyond  the  scope  of  this  work  to  formally  treat  the  effects  of 
advection  on  the  dynamics  of  the  mixed  layer,  and  we  shall  restrict  our  study 
to  steady-state  situations.  However, the  other  two  mechanisms  referred  to 
(penetration  of  radiation  and  planetary  rotation)  will  be  included  in  the  model 
derived  in  the  following  sections. 
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A.  EFFECT  OF  THE  PENETRATING  RADIATION  ON  THE  BUOYANCY  FLUX 

The  simplest  model  for  the  penetration  of  radiation  below  the  sea-surface 
assumes  that  the  net  long-wave  solar  irradiance  (  infra-red  and  red  )  is 
absorbed  at  the  surface,  and  the  short-wave  irradiance  decays  exponentially 
with  depth.  Recent  models  (  Paulson  and  Simpson,  1977  ;  Simpson  and  Dickey, 
1980)  also  allow  for  the  penetration  of  long-wave  radiation.  In  both 
approaches,  the  predicted  vertical  profile  of  absorbed  irradiance  is  sensitive 
to  the  values  of  two  empirical  parameters:  the  vertical  scale  of  penetration, 
which  mainly  depends  on  the  turbidity  of  the  water,  and  the  short-wave 
(blue-green)  fraction  of  the  net  solar  irradiance  incident  on  the  surface. 
Because  there  is  only  a  limited  amount  of  data  on  the  optical  properties  of  the 
oceans,  the  values  of  these  empirical  parameters  will  either  have  to  be 
assumed  or  adjusted  as  tunning  constants.  For  that  reason,  and  also  because 
a  steady-state  model  will  be  applied  to  large  oceanic  regions  having  variable 
optical  properties,  the  simplest  model  will  be  used,  i.,e.,  only  the  short- 
wave radiation  will  be  assumed  to  penetrate  below  the  surface.  The  vertical 
profile  of  the  penetrating  irradiance  will  have  therefore  the  following 
exponential  form 

Q2(z)  =  RQS  ez/X     .  (3.1) 

Here,  the  downward  short-wave  irradiance,  Q2(z),  is  the  radiant  flux  density 

on    a    horizontal    surface    due   to    contributions    from    the    entire    upward 
hemisphere;  Qs  is  the  net  solar  irradiance  at  the  surface,  R  is  the  short-wave 

fraction     of  the  net   solar  irradiance   (    blue-green    ),    z    is   the   vertical 
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coordinate,    positive    upward;    and   A.    is   the    penetration    scale,    assumed 
constant  with  depth. 

Alternatively,  the  long-wave  solar  irradiance,  assumed  to  be  totally 
absorbed  at  the  surface,  is  given  by 

Ql«Q0-RQs,  (3.2) 

where  Q0  =  Qs- Gl^- G^- Qe  is  the  net  solar  irradiance  at  the  surface,  minus 

the  long-wave  back  radiation,  minus  the  sensible  heat  flux,  minus  the  latent 
heat  flux. 

Assuming  horizontal  homogeneity  of  the  variables,  it  can  be  shown  that 
the  effective  heat  flux  Q  contributing  to  the  temperature  flux  over  the  mixed 
layer  will  have  the  form  (  Garwood,  1977  ) 

Q  =  Qr  J  I  dQ2/  dz  -  2/h  J  dQ2/ds  ds  1  dz  (3.3) 

or  Q  =  Q0  -  R  Qs  [  2\/h  -   e  _h/^  (  1  +  2\/h  )].  (3.4) 

For  deep  mixed  layers  (  h/\  >>1)  or,  puting  R  =  0,  expression  (3.4) 
reduces  to  Q  -  Q0,  which  corresponds  to  the  net  heat  flux  being  absorbed  at 

the  surface.  For  shallow  mixed  layers  (  h/\<<! ),  expression  (3.4)  reduces 
to    Q  =  Q0  -  R  Qs.  Under  these  circumstancies  the  effect  of  the  penetrating 

radiation  on  the  effective  heat  flux  Is  maximum. 

It  is  also  apparent  from  expression  (3.4)  that,  in  summertime  conditions 
(  small  h,  large  Qs  )    and  in  clear  oceanic  waters  (  large  A.),    Q    can  be 
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significantly  smaller  than  Q0,  thus  reducing  the  buoyancy  flux  and  Increasing 

the  vertical  extent  of  the  wind  generated  turbulence. 

From  expression  (2.9),  the  effective  vertically  integrated  buoyancy  flux 
will  then  be 

B  =  -aghQ/(  2pCp  ),  (3.5) 

where  Q  is  given  by  (3.4). 

A.  THE  ROTATION  STRESS  MECHANISM 

The  steady-state  turbulent  kinetic  energy  budgets  for  the  three  spatial 
components  (  x  eastward,  y  northward,  z  upward  ),  vertically  integrated 
over  the  mixed  layer  are: 

0  =  Gx  +  Px  +  Qz  h  uV(0)  +  Qy  hTx/p  -  D/3  (3.6a) 

0  =  Gy  +  Py-Qz  h  uV(0)  -D/3  (3.6b) 

0  =  B    +  PZ  -  Qy  hTx/p  -  D/3.  (3.6c) 

Here  Gx  and  Gy  are  the  horizontal  components  of  the  shear  production,  B  is 
the  vertically  integrated  buoyancy  flux,  the  Pj"s  are  the  components  of  the 
vertically  Integrated  turbulent  transport,  and  tx  Is  the  zonal  component  of 
the  surface  wind-stress.  The  Qj's  are  the  two  spatial   components  of  the 

planetary  rotation  and  D  is  the  dissipation,  assumed  to  be  equally  partitioned 
among    components    (local    isotropy    assumed).    Since   the    transport    and 


23 


planetary  rotation  terms  vanish,  the  sum  of  the  three  above  equations  simply 
leads  to  : 

B  +  G-D  =  0.  (2.12) 

An  aditional  process,    not  apparent  in   (2.12),    is   evident  in  equations 

(3.6).  This  is  the  exchange  between  horizontal  and  vertical  turbulent  kinetic 
energy,   in  the  presence  of  the  zonal   wind  stress  tx  and  the   meridional 

component     Qv  of  the  planetary   rotation.    Note   that   the   sign   of   tx   will 

determine  whether  the  term     Qyhtx/p  is   a   source  or   sink   of   horizontal 

turbulent  kinetic  energy.     With  easterly  winds   (  tx  <  0   ),  the  horizontal 

kinetic  energy  is  expected  to  be  enhanced  at  the  expense  of  the  vertical 

kinetic  energy,  with  a  consequent  shallowing  of  the  mixed  layer.  With 
westerly  winds  (  tx  >  0  ),  the  exchange  between  the  vertical  and  horizontal 

components  will  have  the  opposite  sign,  leading  to  an  increase  of  the  vertical 
extent  of  turbulent  mixing.  This  mechanism,  usually  overlooked  in  planetary 
boundary  layer  models,  was  first  examined  by  Garwood  et  aJ.  (1985b),  in  an 
attempt  to  explain  the  unusually  deep  extent  of  mixing  in  the  central  and 
western  Pacific. 

If  the  dissipation  is  parameterized  in  terms  of  the  vertical  average  of  the 
total  turbulent  kinetic  energy,  i.e, 

D=m1<E>3/2,  (3.7) 
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where    mi  and  subsequent  mfs  are  dimensionless  constants,  (2.12)  can  be 
written 

G  +  B-  m\  <  E  >3/2  =  0  .  (3.8) 

Also,    the  vertical    component  Pz    of   the   turbulent   transport    can    be 

parameterized  in  terms  of  the  vertical  average  of  the  TKE.  According  to  Rotta 
(1951 ),  a  first  order  approximation  to  the  transport  term  P2  is 

P2  =  m2<  E  >3/2.  (3.9) 

Substituting  (3.9)  into  (3.6c)  and  combining  with  (3.8),  will  yield  a 
system  of  two  equations  which  can  be  solved  for  the  mixed  layer  depth  h  and 
the  vertical  average  of  the  turbulent  kinetic  energy  <  E  >,  given  the  surface 
boundary  conditions  and  the  meridional  component  of  the  planetary  rotation. 
The  predicted  value  of  h  will  then  reflect  the  interaction  between  the  zonal 
wind  stress  and  the  Earth  rotation. 

C.  EQUILIBRIUM  DEPTH  OF  MIXING 
We  define 


L  =  2pCpu*3/(  agQ  )  (3.10) 

G=m3u*3  (3.11) 

and  <t>  =  QyixCp/(  ag  Q),  (3.12) 
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where  the  dimensionless  parameter  $  is  the  quotient  between  the  conversion 
of  vertical  turbulent  kinetic  energy  by  means  of  the  rotation  mechanism  and 
the  buoyant  damping  of  vertical  turbulent  kinetic  energy.  Then,  a  diagnostic 
equation  for  the  mixed  layer  depth  can  be  derived: 

h  =  C1  L/O+Cz*),  (3.13) 

where  the  non-dimensional  constants  C^  and  C2   are  related  to  the  m{s: 

q  =  m3(m2/mr1/3)/(m2/m1+2/3)  (3.14) 

and  C2  =  2/(m2/m1+2/3).  (3.15) 

Note  that  the  effective  surface  heat  flux  Q  in  equations  (3. 10)  and  (3. 12) 
reflects  the  penetration  of  radiation  below  the  surface,  as  modelled  by 
equation  (3.4).  Thus  the  diagnostic  value  of  h  will  be  dependent  on  the  values 
of  four  tunning  parameters:  the  Cj's,  the  factor  R  and  the  vertical  scale  of 

penetration  X.  Because  Q  is  a  relatively  complex  function  of  h,  a  solution  for 
h  from  (3.4)  and  (3.13)  is  only  possible  by  iterative  techniques.  An 
alternative  method,  which  significantly  reduces  the  computational  effort,  is 
to  substitute  the  Obukhov  length  scale  L  for  h  in  equation  (3.4),  and  then 
solve  for  Q.  This  simplification  is  justified  and  will  be  utilized  in  the  next 
sections,  since  the  primary  goal  of  this  study  is  to  examine  the  gross 
characteristics  of  the  solution  space. 

Substituting  (3.13)  back  into  (3.7),  the  following  expression  for  the 
vertically  integrated  dissipation  is  obtained: 
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D=  mj  <  E  >3/2  =  C2  m3  u*3  ( 1+2*)/l2(1K2*)].  (3.16) 


Physically,  the  dissipation  must  always  be  non-negative,  which  yields 


or 


1+2*  2  0 
♦  2-1/2. 


(3.17) 
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Fig.  3.     The  dlmenslonless  depth  scale  h*  plotted  against  the  rotation  stress 
variable  <D. 
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For    ♦  <  -1/2,   and  in  absence  of  advection,    no  steady-state  situation  is 

possible  and  the  mixed  layer  will  continually  entrain.  Also,  noting  that,  in 

equation    (3.13),    the    denominator    must    be    positive,    the    limiting    case 
expressed  by  (3.17)  suggests  that    C2  <  2    for  all  situations  (since    C2  is 

assumed    to    be    a    constant).     Zeman    and    Tennekes     (1975)     obtained 
n^/m^  1/2,      which      gives      C2  ~  12/7.    Using    the    suggested   value,    a 

nondimensional  mixed  layer  depth  Is  given  by  (  Fig.  3  ) 

h*  =  h/C1L=  1/(  1  +  12*/7  ).  (3.18) 

Note  that  to  the  limiting  value  #  =  -1/2  corresponds  a  theoretical  maximum  of 
h*  =  7. 
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IV.     THE    EQUILIBRIUM    MIXED    LAYER   DEPTH    IN   THE   TROPICAL   ATLANTIC 


A.  THE  BOUNDARY  CONDITIONS 

In  the  following  sections,  the  method  of  application  of  the  steady-state 
mixed  layer  model  to  an  oceanic  macro-scale  domain  is  presented.  The  goal 
of  this  study  is  to  gain  insight  into  the  relative  importance  of  the  physical 
mechanisms  involved,  as  well  as  to  identify  the  regions  and  the  time  of  the 
year  where  the  observed  turbulent  layer  can  be  represented  realistically  by 
such  a  simple  steady-state  model.  This  goal  is  limited  to  some  degree  by  the 
availability  and  suitability  of  experimental  observations  of  the  boundary 
conditions.  Although  a  significant  effort  has  been  made  in  the  last  decade  to 
obtain  accurate  and  detailed  measurements  of  the  surface  boundary 
conditions,  such  data  are  frequently  confined  to  limited  oceanic  areas.  The 
available  climatological  data  for  large  oceanic  areas  are  often  unreliable, 
sparse  and  distributed  irregularly,  depending  on  the  usual  routes  followed  by 
mariners.  Given  these  limitations,  the  tropical  Atlantic  Ocean  is  the  best 
covered  ocean  basin. 

Except  for  the  observed  surface  wind-stress,  where  the  digitized  data  of 
Hellerman  and  Rosenstein  (1983)  were  used,  all  the  other  data  were 
interpolated  from  climatological  atlases.  The  area  chosen  was  the  tropical 
and  equatorial  Atlantic  between  6  degrees  south  and  30  degrees  north. 

1 .  The  Observed  Mixed  Layer  Depth 

The  atlas  of  Robinson  et  at.  (1979)  was  used  to  obtain  the  observed 
monthly  values  of  the  mean  depths  to  the  top  of  the  thermocline  for  the 


29 


equatorial  and  tropical  Atlantic.  The  top  of  the  thermocline  is  therein  defined 
as  the  depth  at  which  the  temperature  Is  1 . 1  degrees  centigrade  less  than  the 
surface  temperature.  Although  this  definition  differs  from  the  usual 
definition  of  "mixed  layer  depth"  (the  vertical  extent  of  turbulent  mixing),  the 
difference  should  not  be  significant  for  the  present  study,  provided  density  is 
mostly  dependent  upon  temperature,  and  the  boundary  layer  Is  well  mixed. 
For  the  cases  where  the  mixed  layer  mean  temperature  decreases 
significantly  with  depth,  the  above  criterion  might  underestimate  the  vertical 
extent  of  mixing. 
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Fig.  4.  Zonal  variation  of  the  mean  mixed  layer  depth  for  the  month  of 
September  (equator).  The  curve  is  the  result  of  an  adjusted  cubic  spline  fit 
over  the  observed  values,  represented  in  the  figure  by  circles. 
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Because  the  contours  displayed  in  the  atlas  have  a  relatively  large 
increment  of  15  meters,  a  cubic  spline  fit  was  made  for  each  parallel  of 
latitude,  and  the  resulting  curve  was  further  adjusted  so  that  the  interpolated 
values  for  each  interval  would  lie  between  extremes.  The  resulting  maximum 
graphical  error  is  estimated  to  be  ±10  meters.  No  indication  is  given  by  the 
authors  about  the  estimated  errors  for  their  analized  depths. 

Fig.  4     shows  the  zonal  equatorial  profile  of  the  mean  mixed  layer 
depth  for  the  month  of  September,  as  computed  by  the  above  technique. 
2.  The  Surface  Heat  Flux 

The  heat  budget  atlas  of  Hastenrath  and  Lamb  (1978)  was  used  to 
provide  best  estimates  of  the  net  surface  heat  flux  and  the  net  surface  solar 
irradiance  for  the  area.  For  the  assumptions  made  and  the  calculation 
details,  the  reader  is  referred  to  the  introductory  part  of  the  atlas. 

The  contour  interval   for  the  surface  heat  flux  is  40  W/m2,   which 

constitutes  a    low  resolution.   In  order  to  minimize  the  resulting  errors,  a 

cubic  spline  surface  was  first  fitted  over  a  grid  of  digitized  values,  and  the 

resulting   interpolated  values  were  further   adjusted  for   each   parallel    of 

latitude.  Thus,  for  each  interpolation  interval,  the  calculated  values  would 

lie  between    extremes.  The  considerable  effort  required  for  this  process  is 

well  justified  because  of  the  sensitivity  of  the  diagnostic  equations  to  small 
deviations  in  the  values  of  Q0,  specially  when  Q0   is  small. 

The  relative  and  absolute  errors  in  the  net  surface  heat  flux,  as 
indicated  by  the  authors,  are  estimated  at  less  than  10  and  more  than 
20  W/m2,  respectively.  The  absolute  error  in  interpolated  values  is 
estimated  to  be  less  than  20  W/m2. 
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As  an  example,  Fig.  7a  shows  the  zonal  equatorial  profile  of  the  mean 
surface  heat  flux  for  the  month  of  September,  as  determined  by  the  above 
technique. 

3.  The  Surface  Wind  Stress 

The  surface  wind  stress  data  are  those  of  Hellerman  and  Rosenstein 
(1983),  for  a  mesh  grid  of  2  degrees  of  latitude  by  2  degrees  of  longitude, 
which  is  suitable  for  our  purpose.  The  two  components  are  calculated  by 
means  of  the  bulk  aerodynamic  formulation: 

tx=  Pa<4)  Ux  (Ux2  +  Uy2  )1/2  (4.1) 

and  Ty=  paC0Uy  (Ux2  +  Uy2  )1/2,  (4.2) 

where  Ux  and  Uy  are  the  two  components  of  the  wind  speed  at  10  meters 
above  the  sea-surface.  The  drag  coefficient,  Cn,   depends  upon  wind  speed 

and  stability,  according  to  the  formulation  of  Bunker  (1976).  The  maximum 
standard  error  in  the  surface  wind-stress  values  for  the  area  considered  is 
estimated  to  be  less  than  0.01  N/m2. 

The  wind  stress  data  are  given  for  each  odd  whole  degree  of  latitude 
and  longitude,  and  in  the  present  scheme  we  will  be  working  with  even 
parallels  of  latitude.  Thus  the  original  data  had  to  be  interpolated.  The 
resulting  interpolation  errors  are  considered  to  be  negligible  for  the  purpose 
of  the  present  study. 

4.  Steadv-State  versus  Time-Dependent  Solutions 

Formally,  a  steady-state  mixed  layer  model  should  only  be  applied  to 
a  real  oceanic  situation  when  entrainment  has  stopped  and  the  top  of  the 
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thermocllne  has  retreated  to  an  equilibrium  level,  dependent  on  the  balance 
between  the  downward  surface  heat  flux  and  the  surface  wind-stress.  On  a 
seasonal  time  scale,  i.e.,  neglecting  daily  variations,  and  for  mid- 
latitudes,  the  above  situation  Is  expected  to  occur  in  spring  and  summer, 
when  the  mean  downward  surface  heat  flux  is  rising  and  the  wind  speed  is 
decreasing.  For  the  equatorial  and  tropical  regions,  however,  this  process 
is  complicated,  since  the  mean  surface  heat  flux  follows  a  complex  annual 
variation,  while  the  wind  speed  has  a  strong  annual  component. 

For  simplicity,  the  diagnostic  mixed  layer  depth  should  be  calculated 
for  the  same  month  throughout  the  chosen  area,  though  the  requirement  of 
steadiness  will  not  be  satisfied  in  some  regions.  A  logical  and  simple  way  to 
solve  this  problem  is  to  plot  the  annual  time  series  of  the  zonal  profiles  of 
the  observed  mixed  layer  depths  for  some  parallels  of  latitude  and  choose 
the  time  of  the  year  which  best  conforms  to  the  condition  of  steadiness  for  all 
latitudes.  Figure  5  shows  these  plots,  as  compiled  and  further  interpolated 
from  Robinsons's  atlas.  The  month  of  September  was  chosen  as  conforming 
best  to  these  conditions.  Except  at  the  equator,  where  the  top  of  the 
thermocline  seems  to  be  still  deepening,  the  other  latitudes  show 
approximate  steady-state  situations,  during  September. 

The  analysis  of  Fig.  5a  reveals  a  strong  annual  component  of  the 
equatorial  mixed  layer  depth,  which  suggests  the  domination  of  the  wind 
stress  in  determining  the  vertical  extent  of  mixing.  In  Fig.  6,  the  annual  time 
series  of  the  net  surface  heat  flux  and  the  wind  speed  are  shown.  Indeed,  a 
qualitative  compairason  of  the  contour  lines  of  both  figures  (5a  and  6) 
reveals  a  much  stronger  correlation  between  the  mixed  layer  depth  and  the 
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Fig.  5.    Annual  zonal  variation  of  the  mean  mixed  layer  depth,  (a)  Equator, 
(b)  5  degrees  north. 
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Fig.  5(cont. ).    Annual  zonal  variation  of  the  mean  mixed  layer  depth 
degrees  north,  (d)  15  degrees  north. 
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Fig.  6.  (a)  Annual  zonal  variation,  in  W/m2,  of  the  net  surface  heat  flux  at 
the  equator,  (b)  Annual  zonal  variation,  in  meters,  of  the  wind  speed  at  the 
equator. 
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wind  speed  than  between  the  mixed  layer  depth  and  the  surface  heat  flux.  In 
particular,  the  wind  speed  maximum  (August  to  October)  and  minimum 
(March  to  April),  coincide  well  with  the  corresponding  relative  extremes  in 
the  mixed  layer  depth.  Similar  agreement  is  not  apparent  between  the  mixed 
layer  depth  and  the  surface  heat  flux,  as  so  far  as  an  equilibrium  state  is 
concerned. 

On  the  other  hand,  one  might  expect  the  annual  variation  of  the  surface 
heat  flux  to  show  a  marked  response  due  to  the  annual  variation  of  the  Sun's 
declination,  i.e.,  to  show  a  strong  biannual  component.  However,  the 
exchange  of  heat  between  the  ocean  and  the  atmosphere  is,  in  the  equatorial 
and  tropical  regions,  significantly  influenced  by  other  variables.  These 
variables,  such  as  cloudiness  and  precipitation,  have  seasonal  variations 
which  are  more  complex  and  difficult  to  explain. 

During  the  months  of  May  to  September, another  Interesting  feature  in 
the  eastern  part  of  the  Equator  is  the  mixed  layer  depth  minimum.  This 
occurs  when  the  surface  heat  flux  is  negative  (upward).  Without  advection, 
an  equilibrium  situation  is  obviously  precluded.  Interestingly,  the  atlas  of 
Robinson  et  at.  (1979)  shows  that  this  is  also  a  zone  of  strong  equatorial 
upwelling  and  surface  divergence.  Thus,  an  advective  mechanism  may 
largely  influence  the  dynamics  of  turbulent  mixing  for  that  region. 

From  the  above  discussion  and  also  from  the  knowledge  that  the 
equatorial  and  tropical  regions  are  zones  of  important  oceanic  heat 
divergence,  it  might  be  concluded  that  horizontal  and  vertical  advection  play 
an  important  role  in  determining  the  vertical  extent  of  mixing  at  the  equator. 
As  lateral  transport  in  the  ocean  diverts  much  of  the  heat  exchanged  with  the 
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atmosphere,  this  effect  would  tend  to  spaclally  homogenize  the  surface  heat 
flux  available  for  the  vertical  turbulent  processes.  Thus,  the  seasonal 
variations  of  the  mean  mixed  layer  depth  would  primarily  reflect  the 
seasonal  variations  of  the  wind  speed. 

B.   METHOD  OF  MODEL  TUNING 

The  models  described  in  sections  II  and  III  shall  be  applied  to  a  set  of 
some  parallels  of  latitude  in  the  tropical  and  equatorial  Atlantic,  from  6 
degrees  south  to  30  degrees  north.  For  each  parallel,  the  spacing  of 
calculated  values  will  be  1  degree  of  longitude.  The  purpose  of  the  following 
procedure  is  to  identify  the  zones  where  the  observed  mixed  layer  is  well 
represented  by  the  above  steady-state  models,  as  well  as  to  analyse  the 
relative  influence  of  the  physical  mechanisms  discussed  in  section  III:  the 
penetration  of  radiation  below  the  surface  and  the  rotation  stress 
mechanism.  For  each  parallel  of  latitude,  a  zonal  profile  of  the  following 
nondimensional  depth  will  be  first  calculated: 

H*  =  H(1+C20)/L  ,  (4.3) 

where  H  is  the  observed  mixed  layer  depth  and  C2  =  12/7.  To  calculate  the 
Obukhov  length  scale  L  and  the  rotation  stress  variable  0,  the  penetration 
parameters  will  be  set  to  ,  respectively,  R=O.Z  and  X  =  6  meters.  Note 
that,  for  each  longitude  point,  H*  is  numerically  equal  to  the  constant  C^.  In 

the  ideal  case  of  a  perfect  agreement  between  the  model's  results  and  the 
observations,  the  resulting  curve  of  H*  for  each  latitude  would  be  a  constant 
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of  order  1.  The  departure  of  the  calculated  values  from  that  ideal  solution 

will  constitute  a  method  of  identifying  the  zones  where  steadiness  is  not 

verified  or  advection  plays  an  important  role.  For  the  situations  where  the 

surface  heat  flux  is  upward,  i.e. ,  where  a  one-dimensional  equilibrium  is  not 
possible,  H*  will  be  set  to  0.  The  described  procedure  also  makes  possible 

the  selection  of  the  areas  where  the  models  will  be  applied.  For  these  areas, 
the  below  defined  diagnostic  mixed  layer  depths  will  then  be  calculated: 

(i)     h0  :  the  equilibrium  mixed  layer  depth,  assuming  R=0,  C2=0.  This 

corresponds  to  the  steady-state  Kraus  and  Turner  diagnostic  depth,  where 

radiation  is  assumed  to  be  totally  absorbed  at  the  surface. 

(ii)      h]  :  the  equilibrium  mixed  layer  depth,  assuming  R=0.2,   X  -  6m, 

C2=0.  This  is  the  same  as  h0,  except  that  radiation  is  assumed  to  penetrate 

below  the  surface,  according  to  the  model  given  by  equation  (3.4). 

(iii)     h£  :  the  equilibrium  mixed  layer  depth,  assuming  R=0.2,  X  =  6m  and 

C2*0.  Here,  the  rotation  stress  mechanism  is  also  included,  according  to  the 
formulation  of  Garwood  eta/.  (1985a). 

The  pertinent  equations  are  given  in  sections  II  and  III,   and  the  following 
constants  are  assumed: 

p  =  1025  Kg/m3 
pa=  1.2Kg/m3 

a  =  2x  10_4C_1 
Cp  =  3890  J  Kg-1  C_1 
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g  =  9.8  m/s2 
0  =7.29x  ICT^s"1. 

For  each  of  the  above  defined  depths  and  for  each  latitude,  the  tunning 
constants  Cj    and  C2  will  be  adjusted,  so  that  the  summation  of  the  squares 

of  the  differences  between  the  diagnostic  mixed  layer  depths  and  the  observed 
mixed  layer  depths  is  a  minimum,  and    C ji  0   and   Oi  C2  iZ.  Following  this 

tuning  process,  a  standard  deviation  will  be  then  computed,  according  to  the 
definition: 

N 

s2  =  I(Hrhi)2/(N-l),  (4.4) 

irl 

where  the  Hj's   are  the   observed   mixed   layer   depths,    the   \\\'s    are   the 

diagnostic  mixed  layer  depths  and  N  is  the  number  of  longitude  points.  The 
standard  deviations  will  serve  as  a  quantitative  verification  of  the  models, 
the  model  assumptions,  and  the  accuracy  of  the  boundary  conditions. 

As  discussed  earlier,  the  rotation-stress  variable  0  is  required 
theoretically  to  be  greater  than  -1/2,  the  value  corresponding  to  a  situation 
of  zero  dissipation.  For  that  reason,  d)  will  be  set  equal  to  -1/2  whenever 
its  calculated  value  is  less  than  -1/2.  During  the  month  of  September,  this  is 
likely  to  occur  in  the  western  Atlantic,  where  the  zonal  surface  wind  stress 
is  large  and  the  surface  heat  flux  is  small. 

The  selection  of  the  factor  R  and  the  penetration  depth  \  is  an 
imprecise  process  because  of  the  paucity  of  data  for  the  optical  properties  of 
the  equatorial  and  tropical  oceans.  Although  the  choice  here  of  values  for  R 
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and  X  was  based  on  a  preliminary  tunning  of  the  model  at  a  specific  site,  it  is 
still  somewhat  arbitrary.  Without  a  more  complete  data  set,  the  inclusion  of 
these  parameters  as  spacially-varying  tuning  parameters  does  not  seem 
justifiable  at  this  time.  Such  a  procedure  would  tend  to  mask  the  effects  of 
other  important  physical  mechanisms,  such  as  horizontal  and  vertical 
advection.  An  alternative  might  be  to  define  two  different  scales  of 
penetration  for  the  long  and  short-wave  radiation,  and  then  to  use  the 
experimental  results  of  Paulson  and  Simpson  (1977).  However,  this  would 
also  complicate  the  problem.  Furthermore,  Paulson  and  Simpson's  data  are 
valid  only  for  a  specific  spatial  and  temporal  situation  in  the  North  Pacific.  In 
the  absence  of  detailed  measurements  of  the  ocean's  optical  properties, 
there  is  no  simple  solution  to  this  problem,  and  the  results  of  the  present 
work  will  be  somewhat  affected. 
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C.   RESULTS 

Appendix  1  shows  the  zonal  variation  of  the  nondimensional  depth  H*  for 

latitudes    0,    10,  20  and  30  degrees  north.  Only  two  particular  cases,   the 
equator  and  10  degrees  north,  will  be  presented  and  discussed  in  the  text. 
1 .   Equator 

Fig.  7  shows  the  equatorial  boundary  conditions  :  the  surface  heat  flux 
Q0,   the  net  surface  solar  irradiance  Qs,  the  total  surface  wind  stress  T,  and 

the  zonal  wind  stress  Tx.  Fig.  8a  shows  the  zonal  variation  of  the 
nondimensional  depth  H#,  computed  using  equation  (4.3).  It  is  apparent  from 

this  result  that  only  a  limited  zone,  between  about    31  degrees  west  and  4 

degrees   west   seems   to   be   in    a    steady-state    balance.    For    some    of   the 

remaining    regions,    the    presence    of   an    upward    total    surface    heat    flux 

precludes  the  possibility  of  an  equilibrium,  causing  the  top  of  the  thermocline 

to   continually  entrain.    We   have    previously   suggested   that   the   turbulent 

processes  near  the  eastern  boundary  are  likely  to  be  dominated  by  equatorial 

upwelling,  which  can  reduce  significantly  the  vertical  extent  of  mixing.  The 
values  from  Fig.  7a  seem  to  confirme  this  hypothesis,  since  the  values  of  H* 

are  less  than  1  near  the  coast.  For  the  western  region,  other  mechanisms 
might  partially  explain  the  departure  of  the  model's  results  from  the 
observations.  The  first  of  those  is  related  to  the  strong  wind-driven  current 
near  the  coast  of  Brazil,  the  "Guiana  current".  This  flow  and  the  associated 
meridional  variations  in  the  surface  heat  flux  might  cause  the  local  dynamics 
to  be  dominated  by  horizontal  advection.  A  second  possibly  contributing 
factor  is  the  influence  of  the  Amazon  River  runoff.  According  to  the  data  on 
the  mean  surface  salinities  in  Robinson's  atlas, this  runoff  extends  to  more 
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Fig.  7.  Mean  surface  boundary  conditions  at  the  equator,  during  September, 
(a)  Net  downward  surface  heat  flux,  (b)  Net  downward  surface  solar 
irradiance. 


43 


|IKN/m2)    > 

12    :            ':/ 

\ 

08/ 
.06 

.04 

.02 

0 

-46        -42       -38        -34        -30        -26        -22        -18        -14        -10        -6 


-2  2  6 

Longitude 


Fig.   7(cont.).    Mean  surface  boundary  conditions   at   the   equator,    during 
September,  (c)   Total  surface  wind  stress,  (d)  Zonal  surface  wind  stress. 
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than  300  nautical  miles  eastward  from  the  coast  and  it  affects  the  upper  few 
tens  of  meters  of  the  water  column.  As  examined  by  Garvine  (1984),  an 
advected  shallow  plume  of  buoyant  water  can  decrease  the  mixed  layer 
depth.  Also,  since  the  water  advected  from  the  Amazon  is  more  turbid  than 
the  average  sea  water,  the  value  chosen  for  the  penetration  scale  A.  is 
probably  too  large  for  the  region,   causing  the  model  to  overestimate  the 

mixed  layer  depth.  Fig.  7a  shows  precisely  this  effect  :  approaching  the 
western  boundary,  H*  falls  markedly  to  about  0.4. 

For  the  equatorial  region  between  31  degrees  west  and  4  degrees  west, 
the  diagnostic  values  were  calculated  for  the  mixed  layer  depths  scales  h0, 

h}  and  h2  .  The  results  were  compared  with  the  observed  mixed  layer  depths. 

Fig.  8  shows  the  zonal  profiles  obtained  for  each  of  the  defined  scales.  The 
vertical  bars   represent  the  effect   of  varying  Q0      by     ±10  W/m^,    which 

indicates  how  sensitive  the  model  is  to  moderate  errors  or  variations  in  the 
heat  flux  data.  The  comparaison  of  Figures  7a  and  8  also  shows  that  h£  is 

significantly  more  sensitive  to  variations  in  Q0  than  are  the  other   depth 

scales,  especially  when  Q0  is  small.  In  determining  a  best  fit  for  h0  and  h^ , 

the  corresponding   standard  deviations   were  s0=   5.4   meters   and  $]  =  6.2 

meters.  These  optimized  values  were  obtained  with  Cp  1 .3  and  Cp  1 .0,  for 
h0andh},  respectively. 

Fig.  9    illustrates  the  tunning  process  of  C^  and  C2,  in  computing  the 
zonal  variation  of  the  depth  scale  h2«  Displayed  here  are  isopleths,  in  m^,  of 

the  sums  of  the  squares  of  the  differences  between  the  diagnostic  and  the 
observed  depths,  as  a  function  of    C]  and   C2    values.  For  each  grid  point  of 

C]  and  C2  values,  a  zonal  profile  of  the  diagnostic  depth  h2  was  calculated, 
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Fig.  8.  Model  results  for  the  equator,  (a)  Zonal  variation  of  the 
nondimensional  mixed  layer  depth  H«.  (b)  Zonal  variation  of  the  Kraus  and 
Turner  diagnostic  mixed  layer  depth  h0.  The  circles  represent  observed 
values,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±10  W/rn2.  The 
shaded  areas  in  (a)  correspond  to  the  zones  where  a  one-dimensional, 
steady-state  balance  is  assumed  impossible. 
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Fig.   8(cont.).   Model   results  for  the  equator,    (c)    Zonal   variation   of   the 
Kraus  and  Turner  diagnostic  mixed  layer  depth  hj.  (d)  Zonal  variation  of  the 

Garwood  diagnostic  mixed  layer  depth  h2-  The  circles  represent  observed 

values,  while  the  small  squares  represent  interpolated  values.  The  vertical 

bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±10  W/m2. 
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as  well  as  the  zonal  summation  of  the  squares  of  the  differences  between 

those  depths  and  the  observed   mixed   layer  depths.    The  arrows  in  Fig. 9 

indicate  the  optimal  values  of  the  constants,   corresponding  to  a  minimum 
standard  deviation  s^  =  3.3  meters. 

There  are  some  interesting  features  to  be  noted  in  the  results: 

(i)  There  is  no  quantitative  improvement  in  the  quality  of  the  results 
for  the  Kraus  and  Turner   model   (h0)   when  the  effect  of  the  penetrating 

radiation  is  included  (hj).  This  is  explained  by  the  fact  that  the  surface  heat 

flux  is  relatively  large  in  that  region,  and  the  vertical  scale  \  of  penetration 

is  much  less  than  the  computed  Obukhov  length  scale. 

(11)  The  zonal  variation  of  the  mixed  layer  depth  near  the  western  limit 
is  better  represented  when  the  rotation  stress  mechanism  is  included  (h^). 

In  that  region,  the  zonal  wind  stress  has  relatively  large  negative  values, 
which  causes  0  to  also  have  relatively  large  negative  values. 

(iii)  The  model  seems  much  more  sensitive  to  the  value  of  Cj  than  to 
the  value  of  C?.  Fig.  9  shows  a  relatively  large  zone,  between  C£= 1 .2  and 
C2=1 .6  where  the  standard  deviation  remains  practically  constant.  This  fact, 
associated  with  the  similarity  of  the  results  between  h0  and  hj ,  might  indicate 

that  the  effect  of  the  rotation  stress  mechanism  is  small, for  the  region  taken 
as  a  whole. 

In  general,  the  revised  model  shows  a  remarkably  good  agreement 
with  the  observations.  However,  and  because  the  Kraus  and  Turner  h0's  are 

already,  per  se,    in  good  agreement  with  the  observed  values,  the  results 
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are  not  conclusive  concerning  the  Importance  of  the  penetrating  radiation  and 

the  rotation  stress  mechanism.  The  eastern  part  of  the  equator, where  the 

observed  zonal  wind  stress  is  positive,  would  constitute  a  suitable  domain  to 

further  test  the  significance  of  the  rotation  stress  mechanism.  There,  the 
diagnostic  depths  h£  would  be  larger  than  the  Kraus  and  Turner  depths  h0  and 

hj.   However,    as  we  have  seen,    other   mechanisms    appear   to   affect   the 
dynamics  of  mixing  In  the  region,  and  will  need  to  be  included  in  such  a  test. 
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Fig.  9.  Variation  of  the  sum  of  the  squares  of  the  differences,  in  m^  .between 
the  observed  equatorial  mixed  layer  depths  H  and  the  diagnostic  mixed  layer 
depths  \\2  with  the  tunning  parameters  C]  and  C2. 
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2.   10  degrees  North 

Fig.    10  shows  the  mean  boundary  conditions  for   10  degress  north 

during  the  month  of  September.  In  Fig.  11a,  the  zonal  variation  of  the 
nondimensional  depth  H*  is  displayed.  This  figure  shows  that  the  zone  which 

seems  to  be  in  a  steady-state  balance  (between  about  38  degrees  west  and  17 

degrees  west)  is  considerably  smaller  than  the  corresponding  zone  for  the 

equator.  Although  the  surface  heat  flux  is  always  positive  for  the  entire 

domain,  the  values  are  usually  smaller  than  at  the  equator.  When  the  effect 

of  the  penetration  of  radiation  is  included,  the  result  is  to  further  reduce  the 

effective  heat  flux  at  the  surface,  as  defined  by  equation  (3.4).  This  effect 

reduces  buoyant  damping  and  precludes  the  possibility  of  an  equilibrium  for  a 

larger  area.  Although  unsteadiness  may  partially  account  for  these  results, 

advection  seems  to  be  a  more  plausible  explanation  for  the  regions  where  the 
effective  heat  flux  is  positive,  since  H*  is  always  less  than  1.  A  meridional 

gradient  is  present  in  the  surface  heat  flux  because  of  the  strong  wind-driven 

current  in  the  western  part  of  the  tropical  Atlantic,   and  this  effect  might 

cause  buoyant  water  to  be  advected  northward,  reducing  the  vertical  extent 
of  mixing.  The  variation  of  H*  near  the  eastern  boundary  is  more  difficult  to 

explain.  From  about  18  degrees  west,  where  it  has  a  value  of  1 ,  H*  continues 

to  increase  eastward,  which  indicates  that  the  observed  depth  becomes 
larger   than   the    diagnostic    depth    h2-    The    analysis    of   the   temperature 

distribution  1n  the  atlas  of  Robinson  et  al.  (1979)  shows  that  this  is  a  region 
of  particularly  strong  coastal  upwelling.  However,  positive  vertical 
advection  can  only  reduce  the  depth   of  the  turbulent   boundary  layer,    as 


50 


50 

\- 

(W/m2) 

3o\ 

/ 

20  \ 

10 

0 

-10 

-58  -54  -50  -46  -42  -38  -34  -30  -26  -22  -18 

Long>tude 


c 

200 

f(W/m2) 

o-o  . . 

190 

180 

or. 

~>»0 

170 

160 

-58  -54         -50         -46         -42         -38         -34  -30         -26         -22  -18 

Longitude 


Fig.  10.  Mean  surface  boundary  conditions  at  10  degrees  north,  during 
September,  (a)  Net  downward  surface  heat  flux,  (b)  Net  downward  surface 
solar  irradiance. 
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Fig.  10(cont).  Mean  surface  boundary  conditions  at  10  degrees  north,  during 
September,  (c)  Total  surface  wind  stress,  (d)  Zonal  surface  wind  stress. 
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stated    by    Garwood    at    al.    (1985b).    Thus,    this    mechanism    cannot    be 

responsible  for  the  difference. 

The  diagnostic  depths  h0,   h^   and  h2  were  calculated  for  the  region 

between  38  degrees  west  and  17  degrees  west.  Fig.    11   shows  their  zonal 

variations.  The  vertical  bars  represent  the  effect  of  varying  Qo  by  ±5  W/m2, 

from  which  it  can  be  concluded  that  the  model  is  much  more  sensitive  here  to 

small  variations  in  the  heat  flux  than  it  was  for  the  equator,  especially  when 

the  penetration  of  radiation  effect  is  included.   This   is   explained  by  the 

already  discussed  fact  than  the  surface  heat  flux  at   10  degrees  north  is 

significantly  smaller  than  it  is  at  the  equator.  In  some  areas,  particularly  in 

the  western  half  of  the  domain,  the  sensitivity  is  so  large  that  the  resulting 

depths  vary  by  a    factor    of   three   when   those    effects    are   considered. 

Considering  that  the   relative   errors   in  the   surface   heat  flux   data,    as 

estimated  by  Hastenrath  and  Lamb  (1978),  are  less  than  10  W/m2,   only  a 

qualitative  interpretation  of  the  results  for  this  latitude  should  be  made. 

For  each  of  the  diagnostic  scales  h0,  h1  and  h2,  the  optimized  tuning 

constants   and  the   corresponding   standard   deviations   have   the   following 
values: 

h0 :  C!  =  2.6,  C2  =  0.0,  s=  11.8m 
h1  :  C1  =  0.8,  C2  =  0.0,  s  =  7.1  m 
h2  :  Cj  =  0.8,  C2  =  0.7,  s  =  5.9  m. 

A   significant   improvement  in  the    quantitative    quality   of   the    results    is 
apparent,  when  the  penetration  effect  is  included  (hi ).  Note  that  the  optimum 
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Fig.  11.  Model  results  for  10  degreees  north,  (a)  Zonal  variation  of  the 
nondimensional  mixed  layer  depth  H«.  (b)  Zonal  variation  of  the  Kraus  and 
Turner  diagnostic  mixed  layer  depth  h0.  The  circles  represent  observed 
values,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±5  W/m2.   The 

shaded  areas  in  (a)  correspond  to  the  zones  where  a  one-dimensional, 
steady-state  balance  is  assumed  impossible. 
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Fig.  1  l(cont).  Model  results  for  10  degrees  north,  (c)  Zonal  variation  of  the 
Kraus  and  Turner  diagnostic  mixed  layer  depth  h^  (d)  Zonal  variation  of  the 
Garwood  diagnostic  mixed  layer  depth  h2.  The  circles  represent  observed 
values,  while  the  small  squares  represent  interpolated  values.  The  vertical 
bars  show  the  effect  of  varying  the  surface  heat  flux  Q0  by  ±5  W/m2. 
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tuned  value  of  C1  was  reduced  from  2.6  to  a   more  plausible  0.8,   which 

approaches  the  value  found  for  the  equator.  In  general,  the  inclusion  of  the 

penetration  effect  tended  to  homogenize  the  diagnostic   mixed  layer  depth 

field.    This    caused    a    better    agreement    between    the    results    and    the 

observations,  particularly  near  the  extremes  of  the  longitudinal  interval.  As 

expected  from  the  analysis  of  the  boundary  conditions  for  the  zonal  wind 

stress,   the  western  part  of  the  domain  is   better    represented  when  the 

rotation   stress    mechanism   is   included.    Also,    a    small    improvement   is 

observed  near  22  degrees  west,  where  the  observed  mixed  layer  depth  has  a 

minimum.  This  is  explained  by  the  zonal  wind  stress  being  positive  in  that 
region.   When  the  eastern  boundary  is  approached,   Tx  becomes  more   and 

more  positive,    causing  the   diagnosed  value   for   h?   to   be    less   than   the 

observed  value. 

Fig.  1 1  illustrates  the  tuning  process  for  h£.  Again,  the  contour  lines, 

in  m^,  of  the  sums  of  the  squares  of  the  differences  between  the  diagnostic 
and  the  observed  depths,   are  displayed  as   a  function  of  C^   and  C;?.   The 

optimal  value  of  C^  =  0.8  has  not  changed  with  the  inclusion  of  the  rotation 
stress  mechanism.  This  tends  to  verify  the  merits  of  the  procedure,  in  the 
sense  that  no  other  physical  mechanisms,  absent  from  the  model's 
formulation,  are  being  compensated  for.  On  the  other  hand,  although  Fig.  11 
shows  that  the  model  is  still  more  sensitive  to  variations  in  Cj  than  to 
variations  in  C2,  the  optimal  value  of  C2  is  much  better  defined  here  than  at 

the  equator.  Also,  this  relative  sensitivity  seems  to  verify  the  importance  of 
the  rotation  effect. 
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Although  a  much  better  representation  of  the  mixed  layer  depth  field 
was  obtained  with  the  revised  model  at  10  degrees  north,  this  verification  Is 
limited  by  the  uncertainty  in  the  boundary  conditions.  As  we  have  already 
noted,  the  model  is  so  sensitive  to  small  variations  in  the  heat  flux  data,  that 
no  detailed  interpretation  of  the  results  is  justifiable. 


1.2    T 


Fig.  12.  Variation  of  the  sum  of  the  squares  of  the  differences, in  m2, between 
the  observed  mixed  layer  depths  H  and  the  diagnostic  mixed  layer  depths  h2 
with  the  tunning  parameters  Cj  and  C2  ( 10  degrees  north). 
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V.    CONCLUSIONS 

Starting    with    the    theoretical    frame    work    of    the    one-dimensional, 

steady-state  Kraus  and  Turner  ( 1967)  model  for  the  surface  ocean  boundary 

layer,  a  revised  theory  was  presented.  This  new  theory  includes  the  effect  of 

the  penetrating  radiation  below  the  surface  on  the  buoyancy  flux,  as  well  as 

the    rotation   stress    mechanism    examined   by   Garwood    ( 1935a, b).    In   the 

presence  of  a  downward  surface  heat  flux  and  a  surface  wind  stress,  the 

penetration    of    radiation    will    tend    to    reduce    the    buoyant    damping    of 

turbulence,    enhancing   vertical    mixing.    The    rotation    stress    mechanism, 

driven  by  the  interaction  between  the  meridional  component    of  the  planetary 
rotation  Qv  and  the  zonal  wind  stress  Tx,   predicts  the  equilibrium   mixed 

layer  depth  to  increase  when  Tx  is  negative  (westward)  and  to  decrease  when 

Tx  is  positive  (eastward). 

To  test  its  response,  and  to  gain  insight  into  the  relative  importance  of  the 
physical  mechanisms  involved,  the  revised  model  was  applied  to  a  set  of 
boundary  conditions  in  the  tropical  Atlantic.  Except  for  the  surface  wind 
stress,  where  digitized  values  were  used,  the  other  boundary  and  oceanic 
conditions  (observed  mixed  layer  depth,  surface  heat  flux  and  surface 
downward  solar  irradiance)  were  interpolated  from  low-resolution 
climatological  atlases.  Since  the  model  is  particularly  sensitive  to  small 
variations  in  the  heat  flux  data,  the  quality  of  the  results  is  considered  to 
have  been  somewhat  affected  by  the  amount  of  uncertainty  in  the  boundary 
conditions.  Also,  in  the  absence  of  detailed  measurements  of  the  optical 
properties  of  the  tropical  oceans,  the  values  of  the  penetration  parameters  X 
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and  R  were  assumed  to  be  constant  over  the  whole  domain.   Given  these 

limitations,  no  detailed  interpretation  of  the  results  is  appropriate  here. 

In  general,  the  response  of  the  revised  model  compares  favorably  with 

the    observations.    The    improvement    of   the    results    when    the    two    new 

mechanisms   are  included  is  very  significant,    particularly  at   10   degrees 

north.  This  was  expected,  in  so  far  as  the  penetration  effect  is  concerned, 

since  at  10  degrees  north  the  mixed  layer  is  relatively  shallower  than  at  the 

equator,  and  the  downward  surface  heat  flux  is  generally  smaller.  For  this 

latitude,  an  important  result  which  seems  to  support  the  need  to  Include  the 
rotation  stress   formulation  is  that  the  optimal  value  of  the  constant   C] 

remains  constant  when  that  process  is  Included  in  the  tuning  of  the  model. 

In  a  recent  paper,  Garwood  <?t  a/.  (1985b)  examined  the  zonal  dependence 
of  the  climatological  mixed  layer  depth  in  the  near-equatorial  Pacific  and 
concluded  that  the  rotation  stress  mechanism  is  a  plausible  explanation  for 
the  deep  mixing  in  the  central  equatorial  Pacific.  Although  that  feature  is  not 
so  pronounced  in  the  tropical  and  near-equatorial  Atlantic,  the  results 
obtained  here  seem  to  support  the  same  conclusion,  especially  in  its  western 
part.  On  the  other  hand,  the  tropical  Atlantic  is  a  relatively  smaller  ocean 
basin,  and  the  influence  of  coastal  physical  processes,  like  coastal 
upwelling,  are  more  likely  to  influence  significantly  the  dynamics  of  the 
boundary  layer.  Indeed,  this  effect  is  apparent  from  the  results,  since  only 
the  central  part  of  the  domain  seems  to  be  in  a  one-dimensional  equilibrium 
state.  Also,  the  sensitivity  of  the  model  to  small  variations  in  the  surface 
heat  flux  shows  that  this  equilibrium  state  is,  as  it  is  in  the  tropical  Pacific, 
a  delicate  one. 
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The  choice  of  the  tropical  Atlantic  as  a  domain  for  application  of  the  model 
was  primarily  motivated  by  the  existence  of  a  suitable  set  of  surface 
boundary  condition  data.  Future  research  on  the  subject  should  extend  the 
domain  to  the  non-tropical  regions  and  should  include  unsteadiness  and 
advection.  Although  the  model  derived  in  this  paper  is  a  simplistic  one- 
-dimensional  representation  of  a  complex  phenomenon,  and  the  quality  of  the 
results  were  probably  affected  by  the  lack  of  detail  of  the  available  data,  it 
is  concluded  here  that  the  physical  mechanisms  of  rotation  stress  and 
penetration  of  radiation  are  important  in  determining  a  steady-state 
equilibrium  depth  of  turbulent  mixing  for  the  tropical  Atlantic. 
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APPENDIX  A  -  MODEL  OUTPUT 
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Fig.   13.  Zonal  variation  of  the  nondimensional  mixed  layer  depth  H*.  Top  : 
equator.   Bottom  :   10  degrees  north.   The  shaded  areas  correspond  to  the 
zones     where     a     one-dimensional,      steady-state 
impossible. 


balance     is     assumed 
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Fig.  14.  Zonal  variation  of  the  nondimensional  mixed  layer  depth  H*.  Top  '.  20 
degrees  north.  Bottom  :  30  degrees  north.  The  shaded  areas  correspond  to 
the  zones  where  a  one-dimensional,  steady-state  balance  is  assumed 
impossible. 
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APPENDIX  B  -  COMPUTER  PROGRAMS 


'  Program  depth  -  calculates  the  steady-state  mixed  layer  depth 

variation  for  a  given  parallel  of  latitude. 
'  Joaquim  Filipe  Gaspar  -  April  1985  -  Naval  Postgraduate  School. 

DEFINT  i-k,n 

DIM  long(200),Qo(200),Qs(200),tau(200),taux(200),depth(200),h(200) 

********************   Variables  ******************************* 

lat  -  latitude  in  degrees  ;  long  -  longitude  in  degrees 

Qo  -  surface  heat  flux  in  W/m?  ;  Qs  -  surface  solar  irradiance  in  W/m^ 

Q  -  effective  heat  flux  in  W/m2 

tau  -  total  wind  stress  in  N/m2  ;  taux  -  zonal  wind  stress  in  N/m^ 

depth  -  observed  mixed  layer  depth  in  meters 

h  -  diagnostic  mixed  layer  depth  in  meters 

L  -  Obukhov  length  scale  in  meters 

Fi  -  rotation  stress  parameter 

R, lb  -  radiation  penetration  parameters 

c  1 , cZ  -  model  tunning  parameters. 

************************************************************** 

Inputdata: 

INPUT  *  Enter  latitude  ';  lat:  lat$=STR$(lat) 
fl$="sh"+lat$  :  OPEN  "1n",*1,f1$ 
f2$="qs"+lat$  :  OPEN  "in*,*2,f2$ 
f3$="tau'+lat$  :  OPEN  "in",#3,f3$ 
f4$=""taux"+lat$  :  OPEN  "1n",*4,f4$ 
f5$="mld'-Hat$   :  OPEN  "in",* 5, f 5$ 

n=0 

WHILE  NOT  EOF(1) 

n=n+l 

INPUT*1,long(n),Qo(n) 

INPUT#2,lg,Qs(n) 

INPUT#3,lg,tau(n) 

INPUT#4,lg,taux(n) 

INPUT#5,lg,depth(n) 
WEND 
CLOSE 
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Inputparameters: 

PRINT  "Longitude  limits  in  files  are  ";  longO  ),long(n) 

INPUT  "Enter  new  longitude  limits  ";  longO,longl 

IF  longO<long(1 )  OR  long1>long(n)  THEN  Inputparameters 

INPUT  "Enter  value  of  CI  ";  cl 

INPUT  "Enter  value  of  C2  ";  c2 
10       INPUT  "Enter  value  of  R  ";  R  :  IF  R<0  OR  R>1  THEN  10 
20       INPUT  "Enter  value  of  lambda  ";  lb  :  IF  lb<=0  THEN  20 

Functions: 

DEF  FNL(tau,Q)=124078!*tau~1.5/Q 

DEF  FNQ(Qo,Qs,L)=Qo-R*Qs*(2*lb/L-(1+2*lb/L)*EXP(-t/lb)) 

DEFFNFi(taux,Q)=144.6*taux*C0S(lat/57.2958)/Q 

Calculate: 
sum=0 
FOR  1=1  TO  n 

IF  long(i)<long0  OR  1ong(i)>long1  THEN  30 
IF  Qo(i)<=0  THEN  h(i)=0  :  c$="*"  :  GOTO  display 
L=FNL(tau(i),Qo(i)) 

Q=FNQ(Qo(i),Qs(i),L):IF  Q<=0  THEN  h(i)=0:c$="*":GOTO  display 
l=FNl(tau(i),Q) 

Fi=FNFi(taux(i),Q)  :  IF  Fi<-.5  THEN  Fi=-.5 
h(i)=d*L/(1+c2*Fi) 
sum=sum+(depth(i)-h(i))"2 
c$="  * 
display: 

PRINT  long(i);c$,h(i) 
30       NEXT  i 

PRINT  "Sum  :  ";  sum 

Choose: 

INPUT  '  Save  results(  1  ),New  flles(Z),  new  parameters(3),  Quit(4)  ";  o 
ON  o  GOTO  Saveresults,  Inputdata,  Inputparameters,  quit 

Saveresults: 

INPUT  "Enter  name  of  output  file  ";  fo$ 

OPEN"out",#5,fo$ 

FORi=1  TOn 

WRITE*5,long(i),h(i) 
NEXT  l 
CLOSE 
GOTO  Choose 

Quit: 
END 
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'  Program  Spline  -  fits  a  cubic  spline  over  a  set  of  (x,y)  points. 

'  Joaquim  Filipe  Gaspar  -  January  1985  -  Naval  Postgraduate  School 

DIM  x(200),y(200),a(200,4),s(200)   :  DEFINT  i,j,n,e 
DIM  a0(200),a1(200),a2(200),a3(200) 

Enterdata: 

INPUT  "Read  disk  file(d)  or  enter  data  from  keyboard(k)  ";  a$ 
IF  a$="k"  THEN  key  ELSE  IF  a$="d"  THEN  disk  ELSE  Enterdata 

disk: 

INPUT  'Enter  name  of  file  to  fit  *;  f$ 

OPEN"in",#l,f$ 

n=1 

WHILE  NOT  EOF(1) 

INPUT#1,x(n),y(n) 

n=n+1 
WEND 
CLOSE*! 
n=n-1  :  GOTO  Endc 

key: 

INPUT  "Enter  number  of  points  ";  n 

FORi=l  TO  n:  PRINT  i,  :  INPUT  x(i),y(i)  :  NEXT  i 

Endc: 

INPUT  "Enter  end  condition  (1,2  or  3)  ";  endcond 

nm2=n-2  :  nm1=n-1 

dxl=x(2)-x(1)  :  dy1  =  (y(2)-y(1))/dxl*6 

defmatrix: 

FORi=1  TO  nm2 

dx2=x(i+2)-x(i+1) 
dy2=(y(i+2)-y(i+1))/dx2*6 
a(1,1)=dxl  :  a(1,2)=2*(dxl+dx2) 
a(i,3)=dx2  :  a(i,4)=dy2-dy1 
dxl=dx2  :  dyl=dy2 
NEXT  i 

ON  endcond  GOTO  solve, end2,end3 

end2: 

a(1,2)=a(l,2)+x(2)-x(l) 
a(nm2,2)=a(nm2,2)+x(n)-x(nml) 
GOTO  solve 
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end3: 

dx1=x(2)-x(1)  :  dx2=x(3)-x(2) 

a(1,2)=(dx1+dx2)*(dx1+2*dx2)/dx2 

a(1,3)=(dx2*dx2-dx1*dx1)/dx2 

dxn2=x(nm1  )-x(nm2)  :  dxnl=x(n)-x(nm1 ) 

a(nm2, 1  )  =  Cdxn2*dxn2-dxn1*dxn1 )/dxn2 

a  ( nm2 , 2 ) = ( dxn  1 +dxn2 ) * ( dxn 1 +2*  dxn2 ) /dxn2 

solve: 

FOR  i=2  TO  nm2 

a(i,2)=a(i,2)-a(i,1)/a(i-1,2)*a(i-1,3) 

a(i,4)=a(i,4)-a(i,1)/a(i-1,2)*a(i-1,4) 
NEXT  i 

a(nm2,4)=a(nm2,4)/a(nm2,2) 
FORi=nm2-1  TO  1  STEP-1 

a(i,4)=(a(i,4)-a(i,3)*a(i+1,4))/a(i,2) 
NEXT  i 
FOR  i=1  TO  nm2 

s(i+1)=a(i,4) 
NEXT  i 

ON  endcond  GOTO  endll ,end22,end33 

end)  1: 

s(1)=0:  s(n)=0:  60TO  coeff 
end22: 

s(1)=s(2)  :  s(n)=s(nm1)  :  GOTO  coeff 
end33: 

s(1)=((dx1+dx2)*s(2)-dx1*s(3))/dx2 

s(n)=((dxn2+dxn1  )*s(nm1  )-dxn1*s(nm2))/dxn2 

coeff: 

FORi=1  TO  nml 

dx1=x(i+1)-x(i) 
a3(i)=(s(i-H)-s(i))/(6*dx1) 
a2(i)=s(i)/2 

aHi)=(y(i+1)-y(i))/dx1-(2*dx1*s(i)+dx1*s(i+1))/6 
a0(i)=y(1) 
NEXTi 

calculate: 

INPUT  "Enter  Interval  betwen  points  ";  dx 
PRINT  "x-limits  on  file  were  :  ",x(  1  ),x(n) 
INPUT  "Enter  new  x-limits  ";  xl  ,x2 
INPUT  "Enter  name  of  output  file  ';  f$ 
OPEN  "out\*2,f$ 
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FOR  i=l  TO  nml 
xx*xl-x(i) 
WHILE  xKx(HI) 

yy=a0(i)+al(i)*xx-»-a2(i)*xx*xx+a3(i)*xx»xx*xx 
WRITER, xl,yy 
xx=xx+dx :  x!=xl+dx 
WEND 
NEXT  i 

xx=x1-x(nml ) 
WHILE  x\<=xZ 

yy=aO(nmt  )+al(nm1  )*xx+a2(i)*xx*xx+a3(i)*xx*xx*xx 
WRlTE*2,xt,yy 
xx=xx+dx :  x1=x1+dx 
WEND 
CLOSE#2 
GOTO  Enterdata 
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'  Program  Graph  -  plots  an  Internally  defined  function  or  a  file  from 

disk. 
'  Jo  a  qui  m  Filipe  Gaspar  -  Jan  1985  -  Naval  Postgraduate  School 

DEFINT  a-z  :  DEFSNGx,y,p:  DIM  xf(700),y(700),yf(700) 
pi=3. 141593  :  YNF=3E+38 


menudef: 

MENU  1,0, 
MENU  1 , 1 
MENU  1,2 
MENU  1,3 
MENU  1,4 

MENU  2,0, 
MENU  2,1 
MENU  2,2 
MENU  2,3 
MENU  2,4 

MENU  3,0, 
MENU  3,1 
MENU  3,2 

MENU  4,0, 
MENU  4,1 

MENU  5,0, 
MENU  5,1 
MENU  5,2 

MENU  6,0, 
MENU  6,1 

MENU  7,0, 
MENU  7,1 
MENU  7,2 
MENU  7,3 
MENU  7,4 

MENU  ON 


Function" 

"Enter  new  function" 
"Use  internal  function" 
"Input  new  file  from  disk" 
"Use  present  file" 
Resolution" 
"Very  High" 
"High" 
"Medium" 
•Low- 
Scaling" 
"Automatic" 
"Enter  scaling  factors" 
X-Limits" 
"Enter  x-limits" 
Grid" 

Yes" 

No" 
Title" 

"Enter  title" 
Run" 

"Start  graph" 
"New  graph" 
"Stop-List" 
"Quit" 


res=4  :  autoscale=1  :  grid=1 

Menuloop: 

IFxrangeTHEN  MENU  7,1,1 
ON  MENU  GOSUB  Menucheck 
GOTO  Menuloop 

Menucheck: 

ON  MENU(0)GOSUB  Function, Resolution, Scaling, Limits, Grid, Title, Start 
RETURN 
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Function: 

ON  MENUU)  GOTO  Newf unction, 01dfunction,Newfi1e,01dfile 

Newfunction: 

MENU  RESET  :  STOP 

01  df  unction: 

foption=0 

MENU  1,2,2:  IF  fflag  THEN  MENU  1,4,1 

RETURN 
Newfile: 

INPUT  "Enter  name  of  file  ";  f$ 

0PEN"in',#1,f$ 

i=0  :  yfmin=YNF  :  yf  max=-YNF 

WHILE  NOT  E0F(1) 

i=i-H  :  INPUT*  1,xf(i),yf(i) 

IF  yf(i)>yfmax  THEN  yfmax=yf(i)  :  imax=i 

IF  yf(i)<yfmin  THEN  yfmin=yf(1)  :  imin=i 

WEND 

CLOSE* 1 

xfmin=xf(  1 )  :  xfmax=xf(i)  :  nx=i 

fflag=1  :  foption=1  :  MENU  1,2,1:  MENU  1,4,2 

RETURN 
Oldfile: 

foption=l  :  MENU  1,2,1  :  MENU  1,4,2 

RETURN 

Resolution: 

res=2~(MENU(1)-1) 

F0Ri=1  T0  4:MENU2,i,1  :  NEXT  i 

MENU  2,MENU(1),2:  RETURN 

Scaling: 

IF  MENU(1)=1  THEN   autoscale=1  :  MENU  3, 1,2  :  MENU  3,2,1  :  RETURN 
autoscale=0:  MENU  3,1,1  :  MENU  3,2,2 
INPUT  "Enter  x-axis  step  ";  xstep 
INPUT  "Enter  minimum  y  ";  ymin 
INPUT  "Enter  maximum  y  ";  ymax 
INPUT  "Enter  y-axis  step  ";  ystep 
IF  foption  THEN  yrange=ymax-ymin 
RETURN 

Limits: 

ON  foption+1  GOTO  limit  1  ,limit2 
limitl: 

INPUT  "Enter  minimum  x  ";xmin  :  INPUT  "Enter  maximum  x  ";  xmax 

xrange=xmax-xmin 

IF  xrange=0  THEN  limitl  ELSE  RETURN 
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limitZ: 

PRIHT  *ln  present  file,  minim  and  max  x's  are  :  *;  xfmin,xfmax 
INPUT  "Enter  minimum  x  ";  xmin  :  IF  xmin<xfmin  THEM  1imit2 
IMPUT  'Enter  maximum  x  ';  xmax  :  IF  xmax>xfmax  THEN  limit2 
xrange=xmax-xmin 
RETURN 

Grid: 

IF  MENU(1)=1  THENgrid=1  :  MENU  5, 1 ,2  :  MENU  5,2, 1  :  RETURN 

grid=0  :  MENU  5,2,2  :  MENU  5,1,1:  RETURN 

Title: 

t$»"  :  INPUT    Enter  title  ";t$  :  RETURN 

Start: 

IF  MENU(  t  )«|  THEN  Graph  :  RETURN 

IF  MENU(  1  )=3  THEN  MENU  RESET  :  STOP 

IF  MENU(  1  )=4  THEN  SAVE  -Graph"  :  SYSTEM 

FOR  i=l  TO  6  :  MENU  i,0, 1  :  NEXT  i  :  MENU  7,2,0 

CLS  :  CALL  TEXTSIZE02)  :  RETURN 


Graph: 

FORi=l  TO  6:  MENU  i, 0,0:  NEXT  i 
CLS 

ONfoption+1  GOTO  graphl  ,graph2 

graphl: 
i=0 

IF  autoscale=1  THEN  ymax=-YNF  :  ymin=YNF 
FOR  graphx=4  TO  488  STEP  res 
i=i+1 

x=xmin+xrange*(graphx-4)/484 
GOSUB  Fx 
IF  autoscale=0THEN30 

IF  y(i)>ymaxAND  ABS(y(1)<>YNF)  THEN  ymax=y(i) 
IF  y(i)<ymin  AND  ABSCy(i)oYNF)  THEN  ymin=y(i) 
30      NEXT  graphx 

yrange=ymax-ymin  :  ni=i-1 
GOTO  setitle 

graph2: 
i=1 
WHILE  xf(i)<xmin 

i=i+1 
WEND 
ii=i :  xmin=xf(i) 
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i=nx 

WHILE  xf(i)>xmax 

WEND 

ifn=i :  xmax=xf(i)  :  xrange=xrnax-xmin 

IF  autoscale  THEN  ymax=yfmax  :  ymin=yfmin  :  yrange=ymax-ymin 

setitle: 

CALL  MOYETOC  100,20)  :  PRINT  t$; 

xlabel: 

CALL  TEXTSIZE(9) 

LINE(4,4)-(4,284)  :  LINE(4,284)-(488,284) 
IF  autoscale  THEN  xstep=xrange/4 
FOR  x=xmin+xstep  TO  xmax  STEP  xstep 

graphx=-484*(xmin-x)/xrange+4 

PSET(graphx,283) 

IF  grid=0THEN40 

FOR  graphy=280  TO  4  STEP  -4 
PSET(graphx, graphy) 

NEXT  graphy 
40  CALL  MOVETOCgraphx-12,295) 

IF  xoxmaxTHEN  PRINT  x; 
NEXTx 

ylabel: 

IF  autoscale  THEN  ystep=yrange/4 
FOR  y=ymin  TO  ymax  STEP  ystep 

graphy=280*(ymin-y)/yrange+284 

PSET(5,graphy) 

IF  grid=0THEN50 

FOR  graphx=9  TO  486  STEP  4 

PSET(graphx, graphy) 
NEXT  graphx 
50  CALL  MOVETO(6,graphy-2) 

PRINT  y; 
NEXT  y 

Plot: 

ON  foption+1  GOTO  plotl , pi ot2 

plotl: 

FOR  1=1  TO  ni  :  IF  (y(i)>ymax  OR  y(i)<ymin)  THEN  NEXT  i 
gy=280*(ymin-y(i))/yrange+284 
gx=4+(i-1  )*res  :  i=0 
FOR  graphx=4  TO  488  STEP  res 
1=1-*- 1 


73 


IF  (y(i)>ymax  OR  y(i)<ymin)  THEN  flagy=0  :  GOTO  5 
f1agy=flagy+1  :  graphy=280*(ymin-y(i))/yrange+284 
IF  flagy>1  THEN  LINE(gx,gy)-(graphx,graphy) 
gx=graphx  :  gy=graphy 
5         NEXT  graphx 

MENU  7,2,1 

RETURN 

plotZ: 

gx=-484*(xrnin-xf(ii))/xrange+4 
gy=280*(ymin-yf(ii))/yrange+284 
FOR  i=ii  TO  ifn 

graphx=-484*(xmin-xf(i))/xrange+4 

graphy=280*(yrnin-yf(i))/yrange+284 

LINE(gx,gy)-(graphx,graphy) 

gx=graphx :  gy=graphy 
NEXTi 
MENU  7,2,1 
RETURN 

Fx: 

y(0=EXP(-x*x) 
RETURN 
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